{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ffd9bbd-3959-46e9-839f-63a386faf005",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c8ddfa93-d8d9-4331-aa6d-8570a190b67f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import struct\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.path import Path\n",
    "import matplotlib as mpl\n",
    "import re\n",
    "import cartopy.crs as ccrs\n",
    "import matplotlib.ticker as mticker\n",
    "\n",
    "\n",
    "%config InlineBackend.figure_format='retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a097f686-d862-4caa-bb1c-5e891d11ba96",
   "metadata": {},
   "outputs": [],
   "source": [
    "## function to sort points counter clockwise (for plotting)\n",
    "def ccw_sort(p):\n",
    "    p = np.array(p)\n",
    "    mean = np.mean(p,axis=0)\n",
    "    d = p-mean\n",
    "    s = np.arctan2(d[:,0], d[:,1])\n",
    "    return p[np.argsort(s),:]\n",
    "\n",
    "## function for great-circle distance between points\n",
    "def great_circle_dist_moon(p1,p2,r=1737.0):\n",
    "    \"\"\"\n",
    "    Calculate great-circle distance between two points\n",
    "    \n",
    "    Input:\n",
    "        p1 = [lon, lat] \n",
    "        p2 = [lon, lat]\n",
    "\n",
    "        r = Radius (Default 1737.0 km) \n",
    "    Output:\n",
    "        d = distance (Default km)\n",
    "    \"\"\"\n",
    "    \n",
    "    lon1 = np.deg2rad(p1[0])\n",
    "    lon2 = np.deg2rad(p2[0])\n",
    "    lat1 = np.deg2rad(p1[1])\n",
    "    lat2 = np.deg2rad(p2[1])\n",
    "    \n",
    "    cent_ang = np.arccos(np.sin(lat1)*np.sin(lat2) + np.cos(lat1)*np.cos(lat2)*np.cos(np.abs(lon1-lon2)))\n",
    "    d = r*cent_ang\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1a6d254-a6f2-4b28-a5df-b98dc033d12a",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Load Data\n",
    "\n",
    "### Thorium\n",
    "https://pds-geosciences.wustl.edu/missions/lunarp/reduced/thoriumhd.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fa783495-b0be-4505-b8c3-a847580eaabd",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/zj/46zkt8ns517_676qkc7f89tm0000gv/T/ipykernel_53065/2391662833.py:2: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support regex separators (separators > 1 char and different from '\\s+' are interpreted as regex); you can avoid this warning by specifying engine='python'.\n",
      "  dat_th = pd.read_csv('../data/thoriumhd.txt', skiprows = 13, sep=',\\s+', lineterminator='\\n', names=[\"lat_min\",\"lat_max\",\"lon_min\",\"lon_max\",\"ppm\"])\n"
     ]
    }
   ],
   "source": [
    "## Load thorium data\n",
    "dat_th = pd.read_csv('../data/thoriumhd.txt', skiprows = 13, sep=',\\s+', lineterminator='\\n', names=[\"lat_min\",\"lat_max\",\"lon_min\",\"lon_max\",\"ppm\"])\n",
    "\n",
    "ppm = dat_th.ppm.values.reshape([360,720])\n",
    "thlon = dat_th.lon_min.values.reshape([360,720])\n",
    "thlat = dat_th.lat_min.values.reshape([360,720])\n",
    "\n",
    "\n",
    "## transform into xarray Dataset\n",
    "thor = xr.Dataset(dict(th_ppm = ([\"x\",\"y\"],ppm)),coords = dict(lat=([\"x\",\"y\"],thlat),lon=([\"x\",\"y\"],thlon))) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d671628-9916-4d72-9dea-93505ceebc32",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7186883f-5fab-4b8c-94a4-89d8cea58ec7",
   "metadata": {
    "tags": []
   },
   "source": [
    "### LROC WAC Mosaic\n",
    "\n",
    "Note: using PDS3 version but recently has been updated to PDS4. Have not yet figure out how to display .TIF file.\n",
    "\n",
    "https://pds.lroc.asu.edu/data/LRO-L-LROC-5-RDR-V1.0/LROLRC_2001/EXTRAS/BROWSE/WAC_GLOBAL/WAC_GLOBAL_E000N0000_004P.TIF\n",
    "\n",
    "First, remove header lines and save to new file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a38542e0-2138-453f-8796-b8d64ca3212e",
   "metadata": {},
   "outputs": [],
   "source": [
    "wacfile = '/work1/moon/WAC/WAC_GLOBAL_E000N0000_004P_nohead.txt'\n",
    "\n",
    "\n",
    "LINES = 720 #rows / lat\n",
    "LINE_SAMPLES = 1440 #collumns /lon\n",
    "SAMPLE_BITS = 32\n",
    "BYTES_PER_PIXEL = int(SAMPLE_BITS/8)\n",
    "\n",
    "llat = np.linspace(-90,90,LINES)\n",
    "llon = np.linspace(-180,180,LINE_SAMPLES)\n",
    "LLON, LLAT = np.meshgrid(llon,llat)\n",
    "\n",
    "pixels = LINES*LINE_SAMPLES\n",
    "fmt = '{endian}{pixels}{fmt}'.format(endian='<',pixels=pixels,fmt='f')\n",
    "#plt.figure(figsize = [20,16])\n",
    "with open(wacfile,'rb') as f:\n",
    "    wacimg = np.array(struct.unpack(fmt,f.read(pixels*BYTES_PER_PIXEL))).reshape(LINES,LINE_SAMPLES)\n",
    "wacimg_corr = np.flip(np.flip(wacimg), axis=1)\n",
    "wacimg_corr_adj = np.roll(wacimg_corr, 1440-210,axis =1) ## manual adjustment for strange offset - this is for display only!!!!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f8d2f0e-7715-49ee-a281-86a0d7d37e7a",
   "metadata": {
    "tags": []
   },
   "source": [
    "### TSTA Results File:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5356c80d-5028-4ff3-a003-da1c0942c38e",
   "metadata": {},
   "outputs": [],
   "source": [
    "full_craterdat = pd.read_csv('../results/tsta_20_200.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4661c6ae-27f1-4291-a93b-20aaaa4e5837",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Constants and Functions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d24b8d46-b9d1-4979-aeaf-7beb0dd70a21",
   "metadata": {},
   "outputs": [],
   "source": [
    "th_c = 0.63 #ppm metzger 1973\n",
    "\n",
    "rho_c = 2550.0 #kg/m**3 Wieczorek 2013  bulk density of highlands crust\n",
    "rho_c_km = rho_c*(1000**3) #kg/km**3\n",
    "\n",
    "D_sc = 18 #km (simple complex transition)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cf6eb7f6-758e-4a77-80be-a39db514980c",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Imbrium (from Neumann 2015)\n",
    "imb_lat = 37.0\n",
    "imb_lon = 341.5 - 360\n",
    "imb_dc = 1321 # main ring \n",
    "dist_from_imb = great_circle_dist_moon([imb_lon,imb_lat],[LLON,LLAT])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "966eab07-047e-4062-b671-f5fc56bf466e",
   "metadata": {},
   "outputs": [],
   "source": [
    "###scenario 1\n",
    "def kreep_over(dat,quad,th_k,X):\n",
    "    \"\"\"\n",
    "    dat = pandas dataset \n",
    "    quad = subset of pandas dataset\n",
    "    th_k = numpy array KREEP thorium concentration\n",
    "    X = numpy array Thickness of KREEP layer\n",
    "    \"\"\"\n",
    "    XX,tt = np.meshgrid(X,th_k)\n",
    "\n",
    "    data = np.zeros((len(th_k),len(X),len(dat.loc[quad])))\n",
    "    error = np.zeros_like(data)\n",
    "    a = 0\n",
    "    for _,rows in dat.loc[quad].iterrows():\n",
    "        x_temp = np.zeros_like(XX) #array height of small paraboloid \n",
    "        vc_temp = np.zeros_like(XX)\n",
    "        vk_temp = np.zeros_like(XX)\n",
    "        x_temp[np.where(XX < rows['h_exc'])] = rows['h_exc'] - XX[np.where(XX < rows['h_exc'])] #array height of small paraboloid \n",
    "        vc_temp = 5*np.pi*rows['r_at']*x_temp**2/2.6 ##volume of small paraboloid\n",
    "        vk_temp = -vc_temp + rows['vols']\n",
    "        data[:,:,a] = (tt*vk_temp + th_c*vc_temp)/rows['vols']\n",
    "        error[:,:,a] = data[:,:,a] - rows['ejecta_m_t']\n",
    "        a += 1\n",
    "\n",
    "\n",
    "    rms_err = np.sqrt(np.mean(error**2,axis = 2))\n",
    "    \n",
    "    return rms_err\n",
    "#### scenario 2\n",
    "def kreep_under(dat,quad,th_k,X):\n",
    "    \"\"\"\n",
    "    dat = pandas dataset \n",
    "    quad = subset of pandas dataset\n",
    "    th_k = numpy array KREEP thorium concentration\n",
    "    X = Depth to KREEP layer\n",
    "    \"\"\"\n",
    "    XX,tt = np.meshgrid(X,th_k)\n",
    "\n",
    "    data = np.zeros((len(th_k),len(X),len(dat.loc[quad])))\n",
    "    error = np.zeros_like(data)\n",
    "    a = 0\n",
    "    for _,rows in dat.loc[quad].iterrows():\n",
    "        x_temp = np.zeros_like(XX) #array height of small paraboloid \n",
    "        vc_temp = np.zeros_like(XX)\n",
    "        vk_temp = np.zeros_like(XX)\n",
    "        x_temp[np.where(XX < rows['h_exc'])] = rows['h_exc'] - XX[np.where(XX < rows['h_exc'])] #array height of small paraboloid \n",
    "        vk_temp = 5*np.pi*rows['r_at']*x_temp**2/2.6 ##volume of small paraboloid\n",
    "        vc_temp = -vc_temp + rows['vols']\n",
    "        data[:,:,a] = (tt*vk_temp + th_c*vc_temp)/rows['vols']\n",
    "        error[:,:,a] = data[:,:,a] - rows['ejecta_m_t']\n",
    "        a += 1\n",
    "\n",
    "\n",
    "    rms_err = np.sqrt(np.mean(error**2,axis = 2))\n",
    "    \n",
    "    return rms_err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98f362b7-4cdb-4410-8390-4ee6111eb4d4",
   "metadata": {},
   "source": [
    "## Analysis of TSTA Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e3ae50fd-2984-461d-a566-b065e4cec7c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "full_craterdat = full_craterdat.loc[full_craterdat['floor_m_t'].notna()& full_craterdat['ejecta_m_t'].notna()] #exclude nan locations\n",
    "\n",
    "## distance from imbrium\n",
    "full_craterdat['dist_from_imbrium'] = great_circle_dist_moon([imb_lon,imb_lat],[full_craterdat['lon'],full_craterdat['lat']])\n",
    "\n",
    "\n",
    "##apparent transient crater\n",
    "full_craterdat['r_at'] = (full_craterdat['rc']**(0.921)*(D_sc/2)**(0.079))/(1.32) #adapted from holsapple 1993 with Rat = Rt/1.3\n",
    "\n",
    "#excavation depth\n",
    "full_craterdat['h_exc'] = 1.3*full_craterdat['r_at']/5.0 #melosh pg. 80\n",
    "\n",
    "full_craterdat['diff'] = full_craterdat['ejecta_m_t'] - full_craterdat['floor_m_t'] \n",
    "full_craterdat['vols'] = np.pi*full_craterdat['h_exc']*full_craterdat['r_at']**2/2\n",
    "\n",
    "\n",
    "\n",
    "full_craterdat = full_craterdat.loc[(full_craterdat['ejecta_m_t'] > 0) & (full_craterdat['floor_m_t'] > 0)]\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "67188e4b-bdd2-4718-a6e9-860cd5223365",
   "metadata": {},
   "outputs": [],
   "source": [
    "###### surface_averages\n",
    "\n",
    "lon_temp = np.arange(-180,185,5)\n",
    "lat_temp = np.arange(-90,91,5)\n",
    "llont,llatt = np.meshgrid(lon_temp,lat_temp)\n",
    "\n",
    "distimbt = great_circle_dist_moon([llont,llatt],[imb_lon,imb_lat])\n",
    "th_k = np.linspace(th_c,30,150) # KREEP thorium concentration\n",
    "\n",
    "\n",
    "\n",
    "llont_r = np.ravel(llont)\n",
    "llatt_r = np.ravel(llatt)\n",
    "imbt_r = np.ravel(distimbt)\n",
    "thor_avg = np.zeros_like(llont_r) + np.nan\n",
    "\n",
    "\n",
    "for i in range(len(llont_r)):\n",
    "    \n",
    "    distarr = great_circle_dist_moon([thlon,thlat], [llont_r[i],llatt_r[i]]) #distance from chosen center point\n",
    "\n",
    "    thor_avg[i] = np.mean(ppm[distarr < 300])\n",
    "    \n",
    "\n",
    "thor_avgs = thor_avg.reshape(*llont.shape) ###roi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6f415a57-065c-4944-a67e-0e36bcf79f2d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/janie/Documents/jupyter/lunar_craters/env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3432: RuntimeWarning: Mean of empty slice.\n",
      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",
      "/Users/janie/Documents/jupyter/lunar_craters/env/lib/python3.9/site-packages/numpy/core/_methods.py:182: RuntimeWarning: invalid value encountered in divide\n",
      "  ret = um.true_divide(\n"
     ]
    }
   ],
   "source": [
    "###thickness of top layer\n",
    "\n",
    "lon_temp = np.arange(-180,185,5)\n",
    "lat_temp = np.arange(-90,91,5)\n",
    "llont,llatt = np.meshgrid(lon_temp,lat_temp)\n",
    "\n",
    "distimbt = great_circle_dist_moon([llont,llatt],[imb_lon,imb_lat])\n",
    "th_k = np.linspace(th_c,30,150) # KREEP thorium concentration\n",
    "\n",
    "\n",
    "\n",
    "llont_r = np.ravel(llont)\n",
    "llatt_r = np.ravel(llatt)\n",
    "imbt_r = np.ravel(distimbt)\n",
    "max_depth = np.zeros_like(llont_r) + np.nan\n",
    "min_thor = np.zeros_like(llont_r) + np.nan\n",
    "\n",
    "for i in range(len(llont_r)):\n",
    "    full_craterdat['dist_from_c'] = great_circle_dist_moon([full_craterdat['lon'],full_craterdat['lat']], [llont_r[i],llatt_r[i]]) #distance from chosen center point\n",
    "\n",
    "    local = full_craterdat['dist_from_c'] < 300\n",
    "    local_r = local & (full_craterdat['diff'] > full_craterdat['back_std'])\n",
    "\n",
    "\n",
    "\n",
    "    max_h = full_craterdat.loc[local_r]['h_exc'].max()\n",
    "\n",
    "\n",
    "    X = np.linspace(0,max_h,1000) # Thickness of KREEP layer\n",
    "    XX,tt = np.meshgrid(X,th_k)\n",
    "\n",
    "    rms_err_r = kreep_over(full_craterdat,local_r,th_k,X)\n",
    "\n",
    "    max_ejr = full_craterdat.loc[local_r]['ejecta_m_t'].max()\n",
    "\n",
    "    max_ej_idx = np.argmin(np.abs(th_k - max_ejr))\n",
    "\n",
    "\n",
    "    rms_min_vals = np.abs(rms_err_r[max_ej_idx,:] - rms_err_r.min()*1.15).ravel()\n",
    "    max_thick_idx = np.argsort(rms_min_vals)[:2].max()\n",
    "\n",
    "    inter = (max_ej_idx,max_thick_idx)\n",
    "    \n",
    "    max_depth[i] = XX[inter]\n",
    "    if np.isnan(XX[inter]):\n",
    "        min_thor[i] = np.nan\n",
    "    else:\n",
    "        min_thor[i] = max_ejr\n",
    "    #d_argsr = np.argmin(rms_err_r, axis = 0)\n",
    "    #t_args = np.arange(0,1000,1)\n",
    "\n",
    "    #XX_r = XX[d_argsr,t_args]\n",
    "    #tt_r = tt[d_argsr,t_args]\n",
    "    #rms_r = rms_err_r[d_argsr,t_args]\n",
    "    ###\n",
    "\n",
    "    ##thor_dist = np.abs(tt_r - max_ejr)\n",
    "    #max_depth[i] = XX_r[np.argmin(thor_dist)]\n",
    "    \n",
    "    #if np.isnan(XX_r[np.argmin(thor_dist)]):\n",
    "    #    min_thor[i] = np.nan\n",
    "    #else:\n",
    "    #    min_thor[i] = tt_r[np.argmin(thor_dist)]\n",
    "depth_arr = max_depth.reshape(*llont.shape)\n",
    "thor_arr = min_thor.reshape(*llont.shape)\n",
    "diffs_r = thor_arr - thor_avgs\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "68438e51-9ea1-42d7-b39d-a422380aaa24",
   "metadata": {},
   "outputs": [],
   "source": [
    "###### buried\n",
    "\n",
    "lon_temp = np.arange(-180,185,5)\n",
    "lat_temp = np.arange(-90,91,5)\n",
    "llont,llatt = np.meshgrid(lon_temp,lat_temp)\n",
    "\n",
    "distimbt = great_circle_dist_moon([llont,llatt],[imb_lon,imb_lat])\n",
    "th_k = np.linspace(th_c,30,150) # KREEP thorium concentration\n",
    "\n",
    "\n",
    "\n",
    "llont_r = np.ravel(llont)\n",
    "llatt_r = np.ravel(llatt)\n",
    "imbt_r = np.ravel(distimbt)\n",
    "max_depth_b = np.zeros_like(llont_r) + np.nan\n",
    "max_thor_b = np.zeros_like(llont_r) + np.nan\n",
    "\n",
    "for i in range(len(llont_r)):\n",
    "    full_craterdat['dist_from_c'] = great_circle_dist_moon([full_craterdat['lon'],full_craterdat['lat']], [llont_r[i],llatt_r[i]]) #distance from chosen center point\n",
    "\n",
    "    local = full_craterdat['dist_from_c'] < 300\n",
    "    local_b = local & (full_craterdat['diff'] < -full_craterdat['back_std'])\n",
    "\n",
    "\n",
    "\n",
    "    max_h = full_craterdat.loc[local_b]['h_exc'].max()\n",
    "\n",
    "\n",
    "    X = np.linspace(0,max_h,1000) # Thickness of KREEP layer\n",
    "    XX,tt = np.meshgrid(X,th_k)\n",
    "\n",
    "    rms_err_b = kreep_under(full_craterdat,local_b,th_k,X)\n",
    "    \n",
    "    max_flb = full_craterdat.loc[local_b]['floor_m_t'].max()\n",
    "    \n",
    "    max_fl_idxb = np.argmin(np.abs(th_k - max_flb))\n",
    "\n",
    "\n",
    "    max_fl_errors = rms_err_b[max_fl_idxb,:] ###rms errors of models with that th_k value\n",
    "    min_dep_idx = np.argmin(max_fl_errors) ###index of minimum error \n",
    "    max_depth_b[i] = X[min_dep_idx]\n",
    "    \n",
    "    if np.isnan(max_depth_b[i]):\n",
    "        max_thor_b[i] = np.nan\n",
    "    else:\n",
    "        max_thor_b[i] = max_flb\n",
    "\n",
    "bdepth_arr = max_depth_b.reshape(*llont.shape)\n",
    "bthor_arr = max_thor_b.reshape(*llont.shape)\n",
    "diffs_b = bthor_arr - thor_avgs\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e28a0f50-1e44-4c32-841a-33d3ed88eae0",
   "metadata": {},
   "outputs": [],
   "source": [
    "km_levels = np.arange(0.5,6.5,0.5)\n",
    "th_levels = np.arange(1.0,12.1,1.0)\n",
    "diff_levels = np.arange(0.0,8.1,1.0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4de5329-f181-4952-9549-77c7b5e8d7c0",
   "metadata": {},
   "source": [
    "# Figure 7"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d498c71-c666-4e75-b0ec-84a32e3da626",
   "metadata": {},
   "source": [
    "##### Find Imbrium Basin Outline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "69d0cfb9-2bae-4140-adf2-c02d16c224bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "## Imbrium (from Neumann 2015)\n",
    "imb_lat = 37.0\n",
    "imb_lon = 341.5 - 360\n",
    "imb_dc = 1321 # main ring \n",
    "dist_from_imb = great_circle_dist_moon([imb_lon,imb_lat],[LLON,LLAT])\n",
    "\n",
    "imb_idx = np.where((dist_from_imb > imb_dc/2 - 4) & (dist_from_imb < imb_dc/2 + 4))\n",
    "\n",
    "imb_idxarr= np.zeros((len(imb_idx[0]),2))\n",
    "\n",
    "\n",
    "imb_idxarr[:,0] = LLON[imb_idx[0],imb_idx[1]]\n",
    "imb_idxarr[:,1] = LLAT[imb_idx[0],imb_idx[1]]\n",
    "\n",
    "imb_sort = ccw_sort(imb_idxarr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb2fd6b4-7f0c-4528-9191-a9bd1b39af9d",
   "metadata": {},
   "source": [
    "##### Find South-Pole Aitken Basin Outline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0f8f73d0-6f2d-4b26-9947-a61da7d8bcf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "spa_lat = -53.0\n",
    "spa_lon = 191.0 - 360.\n",
    "\n",
    "spa_dc = 2400.\n",
    "dist_from_spa = great_circle_dist_moon([spa_lon,spa_lat],[LLON,LLAT])\n",
    "spa_idx = np.where((dist_from_spa > spa_dc/2 - 4) & (dist_from_spa < spa_dc/2 + 4))\n",
    "\n",
    "spa_idxarr= np.zeros((len(spa_idx[0]),2))\n",
    "\n",
    "\n",
    "spa_idxarr[:,0] = LLON[spa_idx[0],spa_idx[1]]\n",
    "spa_idxarr[:,1] = LLAT[spa_idx[0],spa_idx[1]]\n",
    "\n",
    "spa_sort = ccw_sort(spa_idxarr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "03554085-1b95-4024-ae00-3ca4e72b989d",
   "metadata": {},
   "outputs": [],
   "source": [
    "depth_arr_m = np.ma.array(depth_arr,mask= thor_arr < 1.0)\n",
    "bdepth_arr_m = np.ma.array(bdepth_arr,mask= bthor_arr < 1.0)\n",
    "\n",
    "diffs_r_m = np.ma.array(diffs_r,mask= thor_arr < 1.0)\n",
    "diffs_b_m = np.ma.array(diffs_b,mask= bthor_arr < 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "98cbd4cf-6d82-45e1-8260-6c313ce723aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "project = ccrs.Robinson(central_longitude=0, globe=None)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17dc4812-50d2-497f-bbef-0f7d81020201",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure( figsize = [16,15])\n",
    "\n",
    "maps = fig.add_gridspec(nrows=3, ncols=2, left=0.05, right=0.94,wspace = 0.05,hspace = 0.02)\n",
    "cbars = fig.add_gridspec(nrows=3,ncols = 1, left =0.94, right= 1.0, wspace = 0.01, hspace = 0.01)\n",
    "\n",
    "ax_t1 = fig.add_subplot(maps[0,0],projection=project)\n",
    "ax_t1.pcolormesh(LLON,LLAT,wacimg_corr_adj, cmap = 'gray', vmin = -0.03, vmax = 0.2, shading='auto', transform = ccrs.PlateCarree())\n",
    "ax_t1.contourf(llont,llatt,depth_arr_m, levels = km_levels, extend = 'max',cmap = 'turbo',transform = ccrs.PlateCarree())\n",
    "ax_t1.plot(imb_sort[:,0],imb_sort[:,1],\"--\",color = 'white', linewidth = 3, transform = ccrs.PlateCarree())\n",
    "ax_t1.plot(spa_sort[:,0],spa_sort[:,1],\".\",color = 'white', markersize = 2, transform = ccrs.PlateCarree())\n",
    "\n",
    "\n",
    "gl1 = ax_t1.gridlines(linewidth=0.5, color='white', alpha=0.5, linestyle='--', draw_labels=False)\n",
    "gl1.xlocator = mticker.FixedLocator([150, 120, 90, 60,30,0, -30, -60, - 90, -120, - 150])\n",
    "gl1.ylocator = mticker.FixedLocator([-60,-30, 0, 30, 60])\n",
    "\n",
    "ax_t1.set_title('Thickness of Upper Layer', fontsize = 18)\n",
    "\n",
    "ax_t2 = fig.add_subplot(maps[1,0],projection=project)\n",
    "ax_t2.pcolormesh(LLON,LLAT,wacimg_corr_adj, cmap = 'gray', vmin = -0.03, vmax = 0.2, shading='auto', transform = ccrs.PlateCarree())\n",
    "ax_t2.contourf(llont,llatt,thor_arr, levels = th_levels, extend = 'max',cmap = 'viridis',transform = ccrs.PlateCarree())\n",
    "gl2 = ax_t2.gridlines(linewidth=0.5, color='white', alpha=0.5, linestyle='--', draw_labels=False)\n",
    "gl2.xlocator = mticker.FixedLocator([150, 120, 90, 60,30,0, -30, -60, - 90, -120, - 150])\n",
    "gl2.ylocator = mticker.FixedLocator([-60,-30, 0, 30, 60])\n",
    "ax_t2.plot(imb_sort[:,0],imb_sort[:,1],\"--\",color = 'white', linewidth = 3, transform = ccrs.PlateCarree())\n",
    "ax_t2.plot(spa_sort[:,0],spa_sort[:,1],\".\",color = 'white', markersize = 2, transform = ccrs.PlateCarree())\n",
    "\n",
    "ax_t3 = fig.add_subplot(maps[2,0],projection=project)\n",
    "ax_t3.pcolormesh(LLON,LLAT,wacimg_corr_adj, cmap = 'gray', vmin = -0.03, vmax = 0.2, shading='auto', transform = ccrs.PlateCarree())\n",
    "ax_t3.contourf(llont,llatt,diffs_r_m, levels = diff_levels, extend = 'both',cmap = 'YlGnBu_r',transform = ccrs.PlateCarree())\n",
    "gl3 = ax_t3.gridlines(linewidth=0.5, color='white', alpha=0.5, linestyle='--', draw_labels=False)\n",
    "gl3.xlocator = mticker.FixedLocator([150, 120, 90, 60,30,0, -30, -60, - 90, -120, - 150])\n",
    "gl3.ylocator = mticker.FixedLocator([-60,-30, 0, 30, 60])\n",
    "ax_t3.plot(imb_sort[:,0],imb_sort[:,1],\"--\",color = 'white', linewidth = 3, transform = ccrs.PlateCarree())\n",
    "ax_t3.plot(spa_sort[:,0],spa_sort[:,1],\".\",color = 'white', markersize = 2, transform = ccrs.PlateCarree())\n",
    "\n",
    "ax_d1 = fig.add_subplot(maps[0,1],projection=project)\n",
    "ax_d1.pcolormesh(LLON,LLAT,wacimg_corr_adj, cmap = 'gray', vmin = -0.03, vmax = 0.2, shading='auto', transform = ccrs.PlateCarree())\n",
    "cb_km = ax_d1.contourf(llont,llatt,bdepth_arr_m, levels = km_levels, extend = 'max',cmap = 'turbo',transform = ccrs.PlateCarree())\n",
    "ax_d1.set_title('Depth of Buried Layer', fontsize = 18)\n",
    "gld1 = ax_d1.gridlines(linewidth=0.5, color='white', alpha=0.5, linestyle='--', draw_labels=False)\n",
    "gld1.xlocator = mticker.FixedLocator([150, 120, 90, 60,30,0, -30, -60, - 90, -120, - 150])\n",
    "gld1.ylocator = mticker.FixedLocator([-60,-30, 0, 30, 60])\n",
    "ax_d1.plot(imb_sort[:,0],imb_sort[:,1],\"--\",color = 'white', linewidth = 3, transform = ccrs.PlateCarree())\n",
    "ax_d1.plot(spa_sort[:,0],spa_sort[:,1],\".\",color = 'white', markersize = 2, transform = ccrs.PlateCarree())\n",
    "\n",
    "ax_d2 = fig.add_subplot(maps[1,1],projection=project)\n",
    "ax_d2.pcolormesh(LLON,LLAT,wacimg_corr_adj, cmap = 'gray', vmin = -0.03, vmax = 0.2, shading='auto', transform = ccrs.PlateCarree())\n",
    "cb_th = ax_d2.contourf(llont,llatt,bthor_arr, levels = th_levels, extend = 'max',cmap = 'viridis',transform = ccrs.PlateCarree())\n",
    "gld2 = ax_d2.gridlines(linewidth=0.5, color='white', alpha=0.5, linestyle='--', draw_labels=False)\n",
    "gld2.xlocator = mticker.FixedLocator([150, 120, 90, 60,30,0, -30, -60, - 90, -120, - 150])\n",
    "gld2.ylocator = mticker.FixedLocator([-60,-30, 0, 30, 60])\n",
    "ax_d2.plot(imb_sort[:,0],imb_sort[:,1],\"--\",color = 'white', linewidth = 3, transform = ccrs.PlateCarree())\n",
    "ax_d2.plot(spa_sort[:,0],spa_sort[:,1],\".\",color = 'white', markersize = 2, transform = ccrs.PlateCarree())\n",
    "\n",
    "\n",
    "ax_d3 = fig.add_subplot(maps[2,1],projection=project)\n",
    "ax_d3.pcolormesh(LLON,LLAT,wacimg_corr_adj, cmap = 'gray', vmin = -0.03, vmax = 0.2, shading='auto', transform = ccrs.PlateCarree())\n",
    "cb_dif = ax_d3.contourf(llont,llatt,diffs_b_m, levels = diff_levels, extend = 'both',cmap = 'YlGnBu_r',transform = ccrs.PlateCarree())\n",
    "gld3 = ax_d3.gridlines(linewidth=0.5, color='white', alpha=0.5, linestyle='--', draw_labels=False)\n",
    "gld3.xlocator = mticker.FixedLocator([150, 120, 90, 60,30,0, -30, -60, - 90, -120, - 150])\n",
    "gld3.ylocator = mticker.FixedLocator([-60,-30, 0, 30, 60])\n",
    "ax_d3.plot(imb_sort[:,0],imb_sort[:,1],\"--\",color = 'white', linewidth = 3, transform = ccrs.PlateCarree())\n",
    "ax_d3.plot(spa_sort[:,0],spa_sort[:,1],\".\",color = 'white', markersize = 2, transform = ccrs.PlateCarree())\n",
    "\n",
    "cb1 = fig.add_subplot(cbars[0])\n",
    "cb1.axis('off')\n",
    "cbar = plt.colorbar(cb_km,ax = cb1)\n",
    "cbar.set_label('(km)', fontsize=16, color='black')\n",
    "plt.setp(plt.getp(cbar.ax.axes, 'yticklabels'), color='black')\n",
    "cbar.ax.tick_params(labelsize=24)\n",
    "cbar.outline.set_visible(False)\n",
    "\n",
    "\n",
    "\n",
    "cb2 = fig.add_subplot(cbars[1])\n",
    "cb2.axis('off')\n",
    "###min thorium\n",
    "cbar2 = plt.colorbar(cb_th,ax = cb2)\n",
    "cbar2.set_label('Thorium (ppm)', fontsize=16, color='black')\n",
    "plt.setp(plt.getp(cbar2.ax.axes, 'yticklabels'), color='black')\n",
    "cbar2.ax.tick_params(labelsize=24)\n",
    "cbar2.outline.set_visible(False)\n",
    "\n",
    "cb3 = fig.add_subplot(cbars[2])\n",
    "cb3.axis('off')\n",
    "cbar3 = plt.colorbar(cb_dif,ax = cb3)\n",
    "cbar3.set_label('Thorium (ppm)', fontsize=16, color='black')\n",
    "plt.setp(plt.getp(cbar3.ax.axes, 'yticklabels'), color='black')\n",
    "cbar3.ax.tick_params(labelsize=24)\n",
    "cbar3.outline.set_visible(False)\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
